%% 生成一个1Khz信号
    % ft = 1.0;
    Fs = 48000;  % 采样率：1 GHz
    freqs  = [2000 ];  % 混合频率
    amps   = [1.00];
    pn     = 8192  ;
    % ==== 时间轴 ====
    t = 0 : 1/Fs : ((pn-1)/Fs);
    % ==== 合成信号 ====
    sin_wave = zeros(1,pn);
    for i = 1:length(freqs)
        sin_wave = sin_wave + amps(i) * sin(2 * pi * freqs(i) * t );
    end
    M=length(sin_wave);

%% Kaiser窗FIR滤波器设计
    osr = 128;
    fc = 24000;             % 截止频率20kHz
    tband = 24000;           % 过渡带宽50kHz
    beta = 12.26526;           % Kaiser窗形状参数

%% 计算滤波器长度
%     N = round(10*Fs*osr/(2*tband));
    N=80;

%% 设计滤波器
    fir_h   = fir1(N, fc/(Fs*osr/2), 'low', kaiser(N+1, beta));
    rec_h   = fir1(N, fc/(Fs*osr/2), 'low', boxcar(N+1));

%% 分析频率响应
    freqz(fir_h, 1, 2048, Fs*osr);
    title(['Kaiser窗低通FIR: ', num2str(fc/1000), 'kHz截止, N=', num2str(N)]);

%% 过采样插值
    inp_data = zeros(1, length(sin_wave)*osr);  % 插值后总长度
    for i=0 : (length(sin_wave)-1)
        inp_data(1+i*osr) = sin_wave(i+1);
    end

%% 数据滤波
    fir_y=filter(fir_h,1,inp_data);
    rec_y=filter(rec_h,1,inp_data);
    fir_y=fir_y*(osr);
    rec_y=rec_y*(osr);

%% 为了防止FFT计算时频谱泄露,故对数据进行了截取
    rec_yt=rec_y(pn/2+1:pn*osr);
    fir_yt=fir_y(pn/2+1:pn*osr);
    plotCustomSpectrum(sin_wave,pn,Fs,"原始数据");
    plotCustomSpectrum(inp_data,length(inp_data),Fs*osr,"插值数据");
    plotCustomSpectrum(fir_yt,length(fir_yt),Fs*osr,"fir输出数据");
    plotCustomSpectrum(rec_yt,length(rec_yt),Fs*osr,"rec输出数据");


    figure();
    plot(inp_data);
    figure();
    plot(rec_yt);
%%
function [freq, spectrum] = plotCustomSpectrum(data, N,Fs,plot_title)
    % PLOTCUSTOMSPECTRUM 计算信号频谱并绘图
    % 输入：
    %   data      - 输入信号（一维向量）
    %   duration  - 波形时长（秒）
    % 输出：
    %   freq      - 频率向量（Hz）
    %   spectrum  - 频谱幅度

    % 计算采样频率
    % N = length(data);
    % Fs = N / duration; % 采样频率 = 点数 / 时长

    % FFT计算
%     Y = fft(data,N);
%     data = data .*blackman(N);
    Y = fft(data);


    % 幅度谱计算
    P2 = abs(Y / N);          % 双边幅度谱
    P1 = P2(1:floor(N/2)+1); % 单边幅度谱
    P1(2:end-1) = 2 * P1(2:end-1); % 能量校正[3,9](@ref)
    spectrum = P1;

    % 生成频率轴
    freq = Fs * (0:length(spectrum)-1) / N;

    % 转换为dB单位
    spectrum = 20 * log10(spectrum );

    % ==== 绘图 ====
    figure;
    hold on;

    % ==== 主瓣和噪声划分 ====
    mainlobe_mask = freq <= 20000;  % 主瓣定义为 0~20kHz
    noise_mask    = ~mainlobe_mask;

    % 主瓣用红色线绘制
    plot(freq(mainlobe_mask), spectrum(mainlobe_mask), 'r', 'LineWidth', 1.5);
    % 噪声用蓝色线绘制
    plot(freq(noise_mask), spectrum(noise_mask), 'b--', 'LineWidth', 1);
    % 画辅助线标注主瓣界限
    ylimits = ylim;
    line([22000 22000], ylimits, 'Color', [0.5 0.5 0.5], 'LineStyle', ':');

    % 图形设置
    xlabel('频率 (Hz)');
    ylabel('幅度 (dB)');
    title(sprintf('%s\n主瓣 vs 噪声对比（主瓣 ≤ 20KHz）', plot_title));
    legend('主瓣频段 (≤20KHz)', '噪声频段 (>20KHz)', 'Location', 'best');
    xlim([0, min(Fs/2, 100e3)]);
    grid on;
end
%%

% 叠加一个白噪声
%     m=0:M-1;
%     sigma=mean(sin_wave.*sin_wave);%噪声方差
%     sin_wave= sin_wave';
%     v=sqrt(sigma).*randn([M,1]);%噪声系数
%     x=sin_wave+v;%叠加噪声的音频信号